#include "fluid.h"

#include <GL\glew.h>

#define GLFW_NO_GLU
#include <gl\glfw.h>

#include <gl\SOIL.h>

#include <cmath>
#include <omp.h>
#include <basecode\misc_utils.h>
//#include "basecode\shader.h"
#include "distance_function.h"
float roygbiv[256][3];


	FluidSolver::FluidSolver()
	{
		N = 1024*2;
		G             = (real) 0.005;				// Gravitational Constant for our simulation
		spacing       = 3.0;						// Spacing of particles
		k			  = (real) (spacing / 500.0);	// Far pressure weight
		k_near		  = k*10;						// Near pressure weight
		rest_density  = 10.5;							// Rest Density
		sigma         = -0.005;							// visco 1
		beta          = -0.004;							// visco 2
		r = (real) (spacing*4.5);					// Radius of Support
		rsq = r*r;									// ... squared for performance stuff
		inv_r = 1.0f / r;
		SIM_W         = 600;						// The size of the world
		paused = false;

		attracting = 0;
		repelling  = 0;
		attractor.x = -SIM_W;
		attractor.y = -SIM_W;
		do_viscosity_step = false;

		grid_width = (int)ceil( SIM_W * inv_r );
		grid_height = (int)ceil( SIM_W * inv_r );
		grid = std::vector< std::vector<particle*> > (grid_width,grid_height); 
		
		initParticles();
		generateColorLUT();

//#define USE_VBO
#ifdef USE_VBO
		glGenBuffers(1, &VertexVBOID);
		glBindBuffer(GL_ARRAY_BUFFER, VertexVBOID);
		glBufferData(GL_ARRAY_BUFFER, verts.size()*sizeof(glm::vec3), NULL, GL_DYNAMIC_DRAW); // GL_STREAM_DRAW - per render update, GL_DYNAMIC_DRAW per frame
		glBindBuffer(GL_ARRAY_BUFFER, 0);
#endif
		//colorShader = new Shader();
		//colorShader->addStage(std::pair<GLenum, const GLchar*>(GL_VERTEX_SHADER, "color.vert"));
		//colorShader->addStage(std::pair<GLenum, const GLchar*>(GL_FRAGMENT_SHADER, "color.frag"));
		//colorShader->install();
		//colorShader->end();

		int SIZ = (int)SIM_W;
		distance_field = std::vector< std::vector<float> > (SIZ, SIZ); 
		std::vector< std::vector<bool> > map = std::vector< std::vector<bool> > (SIZ, SIZ); 
		std::vector< std::vector<glm::vec2> > borderpoints = std::vector< std::vector<glm::vec2> > (SIZ, SIZ); 
		int image_width, image_height, channels;
		unsigned char *image_data = SOIL_load_image( "hit.png", &image_width, &image_height, &channels, SOIL_LOAD_RGB);
		if( image_data )
		{
			for(int i=0; i<SIZ; i++)
			{
				for(int j=0; j<SIZ; j++)
				{
					int idx = (j*image_width+i)*3;
					if( image_data[idx + 0] == 0 )
					{
						map[i][j] = true;
						attractors.push_back( glm::vec2(i,SIM_W-j) );
					}else{
						map[i][j] = false;
					}
				}
			}

		/*	glEnable(GL_TEXTURE_2D);
			glGenTextures(1, &distance_field_tex);
			glPixelStorei (GL_UNPACK_ALIGNMENT, 1);
			glTexParameteri( GL_TEXTURE_2D,GL_TEXTURE_MIN_FILTER, GL_LINEAR );
			glTexParameteri( GL_TEXTURE_2D,GL_TEXTURE_MAG_FILTER, GL_LINEAR );
			glTexImage2D(GL_TEXTURE_2D, 0, 3, image_width, image_height, 0,
				GL_RGB, GL_UNSIGNED_BYTE, image_data);

		*/	
			SOIL_free_image_data( image_data );
		}

		//distance_field_tex = SOIL_load_OGL_texture("hit.bmp", GL_RGB, distance_field_tex, SOIL_LOAD_AUTO );
		DistanceFunction::dead_reckoning( map, distance_field, borderpoints, 0, 0, SIZ, SIZ);
		
	}


	void FluidSolver::render()
	{
		real xs = 1;
		real ys = SIM_W+1;

		static float _minVelocity = 0;
		static float _maxVelocity = 0;

		float mindens = -rest_density;

		glColor3f(1,1,1);
		glDisable(GL_POINT_SPRITE);
		glDisable(GL_BLEND);
		glDisable(GL_TEXTURE_2D);
		glPointSize(1);
		glDepthMask(GL_TRUE);
//#define DRAWFIELD
#ifdef DRAWFIELD
		glBegin(GL_POINTS);
		for(int i=0; i<SIM_W; i++)
		{
			for(int j=0; j<SIM_W; j++)
			{
				float c = 1.0f - distance_field[i][j];
				glColor3f(c,c,c);
				glVertex2f( i, j);
			}
		}
		glEnd();
#endif
		//glEnable(GL_TEXTURE_2D);	
		//glDisable(GL_BLEND);	
		//glBindTexture(GL_TEXTURE_2D, distance_field_tex );

		//float sx = SIM_W;
		//float sy = SIM_W;
		//glColor3f(1,1,1);
		//glBegin (GL_QUADS);
		//glTexCoord2f (0.0f, 0.0f); glVertex2f (0.0f, 0.0f);
		//glTexCoord2f (1.0f, 0.0f); glVertex2f (sx, 0.0f);
		//glTexCoord2f (1.0f, 1.0f); glVertex2f (sx, sy);
		//glTexCoord2f (0.0f, 1.0f); glVertex2f (0.0f, sy);
		//glEnd ();

		//glDisable(GL_TEXTURE_2D);

		//glEnable(GL_POINT_SMOOTH);
		//glEnable(GL_POINT_SPRITE);
		//glEnable(GL_TEXTURE_2D);	
		//glPointSize(r*2);
		//glEnable( GL_BLEND);
		//glBlendFunc (GL_SRC_ALPHA, GL_ONE_MINUS_CONSTANT_COLOR); // used with glBlendColorEXT, asserts default state is glBlendColor(0,0,0,1);
		//glTexEnvf(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_MODULATE);
		//glDepthMask(GL_FALSE);


		if( particlesCreated > 0 )
		{
			double before_upload = glfwGetTime();
#ifdef USE_VBO
			verts.clear();

			for(int i=0; i<particlesCreated; i++){
				glm::vec3 p = glm::vec3( particles[i].pos.x, SIM_W - particles[i].pos.y, particles[i].rho  );
				verts.push_back( p );
			}

			
			//glEnableClientState(GL_VERTEX_ARRAY);
			//glVertexPointer( 2, GL_FLOAT, 0, &verts[0] );
			//glDrawArrays(GL_POINTS, 0, verts.size() );
			//glDisableClientState(GL_VERTEX_ARRAY);

			glBindBuffer( GL_ARRAY_BUFFER, VertexVBOID );
			glEnableClientState(GL_VERTEX_ARRAY);
			#define BUFFER_OFFSET(i) ((char *)NULL + (i))
			glVertexPointer(2, GL_FLOAT, sizeof(glm::vec3), BUFFER_OFFSET(0));   //The starting point of the VBO, for the vertices

			//int attrib_loc = glGetAttribLocation(colorShader->getProgram(), "vVtxPosition");
			//glVertexAttribPointer(attrib_loc, 3, GL_FLOAT, GL_FALSE, sizeof(glm::vec3), &verts[0].x );
			//glEnableVertexAttribArray(weightPosition);

			glBufferData(GL_ARRAY_BUFFER, verts.size()*sizeof(glm::vec3), &verts[0].x, GL_DYNAMIC_DRAW); // Use insead of sub, cause num particles changes...
			//glBufferSubData(GL_ARRAY_BUFFER, 0, verts.size()*sizeof(glm::vec2), &verts[0].x);
			glDrawArrays(GL_POINTS, 0, verts.size() );
			glDisableClientState(GL_VERTEX_ARRAY);
#else
					glBegin(GL_POINTS);
					for(int i=0; i < particlesCreated; ++i)
					{
			//#define velocitycolor
			#ifdef velocitycolor
						float velMagnitude = glm::length( particles[i].vel );
						if( velMagnitude < _minVelocity )
						{
							_minVelocity = velMagnitude;
						}
						if( velMagnitude > _maxVelocity )
						{
							_maxVelocity = velMagnitude * 0.86f;
						}
						float finalC = ( (_minVelocity+velMagnitude)/_maxVelocity  );
						int c = std::min( std::max( (int) (170.0f - finalC*170.0f), 0 ), 170 );
						glColor3fv( roygbiv[c] ); // same as glColor3f( roygbiv[c][0], roygbiv[c][1], roygbiv[c][2] );	
			#else  // Use pressure for color
						float opacity = ( particles[i].rho ); // try log10()
						//float f = 2.0f * particles[i].force.len();
						//opacity-=mindens;
						opacity/=(rest_density-mindens);
						float finalC = opacity*3;
						//glColor4f( f, f, opacity, opacity);
						int c = std::min( std::max( (int) (170.0f - finalC*170.0f), 0 ), 170 );
						//glColor3fv( roygbiv[c] ); // same as glColor3f( roygbiv[c][0], roygbiv[c][1], roygbiv[c][2] );		
						if( particles[i].color < 0.1f )
						{
							glColor3f(1,0,0);
						}else{
							glColor3f(0,0,1);
						}
			#endif
						glVertex2f( (float) (xs + particles[i].pos.x) , (float) (ys + -particles[i].pos.y) );
					}
					// Done drawing points
					glEnd();
#endif

			uploadTime = glfwGetTime() - before_upload;

			// To disable it all
			glColor3f(1,1,1);
			glDisable(GL_POINT_SPRITE);
			glDisable(GL_BLEND);
			glDisable(GL_TEXTURE_2D);
			glPointSize(1);
			glDepthMask(GL_TRUE);
			glBegin(GL_POINTS);
			for(size_t i=0; i<attractors.size(); i++)
			{
				glVertex2f( attractors[i].x, SIM_W - attractors[i].y );
			}
			glEnd();

			int gx = (int)attractor.x;
			int gy =  (int) (SIM_W -attractor.y);
			if( gx > 0 && gy > 0 && gx < (int)(SIM_W-1) && gy < (int)(SIM_W-1) )
			{

				glm::vec2 gradient = glm::vec2( distance_field[gx-1][gy] - distance_field[gx+1][gy], distance_field[gx][gy+1] - distance_field[gx][gy-1] );
				glm::vec2 finale = attractor + gradient * 100.0f;
				glBegin(GL_LINES);
				glColor3f(1,0,0);
				glVertex2f( attractor.x, SIM_W - attractor.y );
				glColor3f(0,1,0);
				glVertex2f( finale.x, SIM_W - finale.y );
				glEnd();
			}



		}
	}


	void FluidSolver::update()
	{
		if(!paused)
		{
			advance_particles();
			update_grid();
			calc_density();
			calc_pressure();
			if (do_viscosity_step) calc_viscosity();
		}

}

void FluidSolver::advance_particles(){
	// UPDATE
	//
	// This modified verlet integrator has dt = 1 and calculates the velocity
	// For later use in the simulation. 

	// For each particles i ...	
	for(int i=0; i < particlesCreated; ++i)
	{			
		// Apply the currently accumulated forces
		particles[i].pos += particles[i].force;

		// Restart the forces with gravity only. We'll add the rest later.
		particles[i].force = glm::vec2(0, -G);

		// Calculate the velocity for later.
		particles[i].vel = particles[i].pos - particles[i].pos_old;

		// (damping has been removed, since there can never be instability issues
		//  in a verlet based physics simulation, and since it was an ugly hack :-)

		// Normal verlet stuff
		particles[i].pos_old = particles[i].pos;
		particles[i].pos += particles[i].vel;

		// If the particle is outside the bounds of the world, then
		// Make a little spring force to push it back in.
		if(particles[i].pos.x < r) particles[i].force.x -= (particles[i].pos.x-r ) / 2;
		if(particles[i].pos.x >  SIM_W-r) particles[i].force.x -= (particles[i].pos.x - (SIM_W-r)) / 2;
		if(particles[i].pos.y < (0+r) ) particles[i].force.y -= (particles[i].pos.y - (0+r) ) / 2;
		if(particles[i].pos.y > SIM_W-r) particles[i].force.y -= (particles[i].pos.y - SIM_W+r) / 2;


		// Handle mouse interaction. 

		glm::vec2 p_to_attr =  (particles[i].pos - attractor);
		real attr_dist2 = glm::dot(p_to_attr, p_to_attr); // Dot with self is the same as dist squared
		real attr_l = 32;

		//  make mouse attract particles
		if( attracting )
			if( attr_dist2 < attr_l*attr_l && attr_dist2 > 16 ){
				glm::normalize(p_to_attr);
				particles[i].force -= p_to_attr * 0.01f;// * (1.0f / 64.0f );
			}

//#define USE_ATTRACTOR_PARTICLES
#ifdef USE_ATTRACTOR_PARTICLES
		for(int j=0; j<attractors.size(); j++ )
		{
			glm::vec2 p_to_attr =  (particles[i].pos - attractors[j]);
			real attr_dist2 = glm::dot(p_to_attr, p_to_attr); // Dot with self is the same as dist squared
			real attr_l = 2.5f*2.5f;

			if( attr_dist2 < attr_l*attr_l )
			{
				glm::normalize(p_to_attr);
				particles[i].force += p_to_attr * 0.01f;// * (1.0f / 64.0f );
			}
		}
#else
			int gx = (int) (particles[i].pos.x);
			int gy = (int) (SIM_W - particles[i].pos.y);

			glm::vec2 gradient = glm::vec2( distance_field[gx-1][gy] - distance_field[gx+1][gy], distance_field[gx][gy+1] - distance_field[gx][gy-1] );			
			float grad_strength_square = glm::dot(gradient,gradient);
			
			if( grad_strength_square > 0.0f )
			{
				particles[i].force -= 4.0f * gradient;
			}
			
#endif

			//  make mouse repel particles
			if( repelling )
				if( attr_dist2 < attr_l*attr_l ){
					particles[i].force += (particles[i].pos - attractor) * (1.0f / 64.0f);
				}

			// Reset the nessecary items.
			particles[i].rho = 0;
			particles[i].rho_near = 0;
	}
}


void FluidSolver::update_grid()
{
	for(int x=0; x<grid_width; x++){
		for(int y=0; y<grid_height; y++){
			grid[x][y] = NULL;
		}
	}

	for(int i=0; i < particlesCreated; ++i){
		particle *p = &particles[i];
		int gridx = (int) floor( p->pos.x * inv_r );
		int gridy = (int) floor( p->pos.y * inv_r );
		particle *successor = grid[gridx][gridy];
		grid[gridx][gridy] = p;
		p->succ = successor;
	}
}

void FluidSolver::calc_density()
{
	// Calculate the density by basically making a weighted sum
	// of the distances of neighboring particles within the radius of support (r)
	int dirx[] = { -1, -1, +1, +1, -1, +1, 0, 0, 0 };
	int diry[] = { -1, +1, -1, +1, 0, 0, -1, +1, 0 };

	for(int i=0; i < particlesCreated; ++i)
	{
		particle *p1 = &particles[i];
		int int_x = (int) floor(particles[i].pos.x * inv_r);
		int int_y = (int) floor(particles[i].pos.y * inv_r);

		for(int j=0; j<9; j++)
		{
			int gridx = int_x + dirx[j];
			int gridy = int_y + diry[j];

			if( gridx < 0 || gridy < 0 || gridx >= grid_width || gridy >= grid_height ) continue;

			particle *p2 = grid[gridx][gridy];

			while( p2 != NULL )
			{
#pragma region CALC_DENSITY
				if( p1->id > p2->id ) // dont check self or same twice
				{
					// The vector seperating the two particles
					glm::vec2 rij = p2->pos - p1->pos;

					// Along with the squared distance between
					real rij_magnitudesquared = glm::dot( rij, rij );         

					// If they're within the radius of support ...
					if(rij_magnitudesquared < rsq)
					{
						// Get the actual distance from the squared distance.
						real rij_len = sqrtf(rij_magnitudesquared);

						// And calculated the weighted distance values
						real q = (real) 1.0f - rij_len / r;
						real q2 = q*q;
						real q3 = q2*q;

						p2->rho += q2;
						p2->rho_near += q3;
						particles[i].rho      += q2;
						particles[i].rho_near += q3;

					} // eof in range
				} // eof is >
#pragma endregion CALC_DENSITY
				p2 = p2->succ;
			} // while!null
		} // for grid

	} // for each particle...
}


void FluidSolver::calc_pressure()
{
	// PRESSURE
	// Make the simple pressure calculation from the equation of state.
	#pragma omp parallel
	for(int i=0; i < particlesCreated; ++i)
	{
		particles[i].press = k * (particles[i].rho - rest_density);
		particles[i].press_near = k_near * particles[i].rho_near;
	}

	int dirx[] = { -1, -1, +1, +1, -1, +1, 0, 0, 0 };
	int diry[] = { -1, +1, -1, +1, 0, 0, -1, +1, 0 };

	// PRESSURE FORCE
	// We will force particles in or out from their neighbors
	// based on their difference from the rest density.
	for(int i=0; i < particlesCreated; ++i)
	{
		particle *p1 = &particles[i];
		int int_x = (int) floor(particles[i].pos.x * inv_r);
		int int_y = (int) floor(particles[i].pos.y * inv_r);

		// For each of the neighbors
		for(int j=0; j<9; j++)
		{
			int gridx = int_x + dirx[j];
			int gridy = int_y + diry[j];

			if( gridx < 0 || gridy < 0 || gridx >= grid_width || gridy >= grid_height ) continue;

			particle *p2 = grid[gridx][gridy];
			while( p2 != NULL )
			{
				if( p1->id > p2->id ) // only compare particles once. also makes sure we dont compare same. of course.
				{
					glm::vec2 rij = p2->pos - p1->pos;
					real rij_magnitudesquared = glm::dot( rij, rij );         
					
					if(rij_magnitudesquared < rsq)
					{
						real rij_len = sqrtf(rij_magnitudesquared);

						// And calculated the weighted distance values
						real q = (real) 1.0 - rij_len / r;
						real q2 = q*q;

						// calculate the force from the pressures calculated above
						real dm = (p1->press + p2->press) * q +
							(p1->press_near + p2->press_near) * q2;

						// Get the direction of the force
						//vec2f D = rij.normal() * dm;

						real forceMag = (real) 1.0/rij_len * dm;
						glm::vec2 D = rij * forceMag;    //  OPTIMISATION!! here we save a sqrt() call.
						//dX += D;
						p2->force += D;
						p1->force -= D;
					} // eof in range
				}
				p2 = p2->succ;
			}
		} // eof grid loop
		//p1->force -= dX;
	} // eof force-loop
}
void FluidSolver::calc_viscosity()
{
	int dirx[] = { -1, -1, +1, +1, -1, +1, 0, 0, 0 };
	int diry[] = { -1, +1, -1, +1, 0, 0, -1, +1, 0 };

   // VISCOSITY
	for(int i=0; i < particlesCreated; ++i)
	{
		particle *p1 = &particles[i];
		int int_x = (int) floor(particles[i].pos.x * inv_r);
		int int_y = (int) floor(particles[i].pos.y * inv_r);

		for(int j=0; j<9; j++)
		{
			int gridx = int_x + dirx[j];
			int gridy = int_y + diry[j];

			if( gridx < 0 || gridy < 0 || gridx >= grid_width || gridy >= grid_height ) continue;

			particle *p2 = grid[gridx][gridy];
			while( p2 != NULL )
			{         
				if( p1->id > p2->id ) // only compare particles once. also makes sure we dont compare same. of course.
				{
					real rij_x = p2->pos.x - p1->pos.x;
					real rij_y = p2->pos.y - p1->pos.y;
					real len = sqrt( rij_x*rij_x + rij_y*rij_y );
					real q = len / r;
	     
					 real rijn_x = (rij_x / len);
					 real rijn_y = (rij_y / len);
					 // Get the projection of the velocities onto the vector between them.
					 real vel_x = (p1->vel.x - p2->vel.x);
					 real vel_y = (p1->vel.y - p2->vel.y);
					 //double u = dpr(rijn, vel);
					 real u = rijn_x * vel_x + rijn_y * vel_y;
	         
					 if(u > 0.0f)
					 {
						// Calculate the viscosity impulse between the two particles
						// based on the quadratic function of projected length.
						real magnitude = (1.0f - q) * (p2->sigma * u + p2->beta * u*u);
						real I_x = magnitude * rijn_x;
						real I_y = magnitude * rijn_y;
			       
						// Apply the impulses on the two particles
						/*
						particles[ni_idx].vel.x -= I_x * 0.5;
						particles[ni_idx].vel.y -= I_y * 0.5;
			            
						particles[nj_idx].vel.x += I_x * 0.5;
						particles[nj_idx].vel.y += I_y * 0.5;
						*/
			            
						// Using verlet, so move positions forward instead with force.
						p1->force.x -= I_x * 0.5f;
						p1->force.y -= I_y * 0.5f;
			            
						p2->force.x += I_x * 0.5f;
						p2->force.y += I_y * 0.5f;
					 }
				}

				p2 = p2->succ;
			} // while( p2 != NULL )
		}
	}
}



void FluidSolver::activateAttractor( int mx, int my )
{
	attracting = 1;
	attractor = glm::vec2( (real )mx, -my + SIM_W);	
}

void FluidSolver::activateRepelling( int mx, int my )
{
	repelling = 1;
	attractor = glm::vec2( (real )mx, -my + SIM_W);	

	/*int gx = (int)attractor.x;
	int gy = SIM_W - (int)attractor.y;
	if( gx > 0 && gy > 0 && gx < SIM_W-1 && gy < SIM_W-1)
	{
		glm::vec2 gradient = glm::vec2( distance_field[gx-1][gy] - distance_field[gx+1][gy], distance_field[gx][gy+1] - distance_field[gx][gy-1] );			
		float value = distance_field[gx][gy];
		float grad_strength = glm::length( gradient );
		printf("gxgy=%ix%i strength = %f, value:%f \n" ,gx,gy, grad_strength, value );
	}*/

}

void FluidSolver::deactivateAttractor()
{
	attracting = 0;
	repelling = 0;
}

void FluidSolver::addParticles( int mx, int my )
{
	attractor.x = (real) mx;
	attractor.y = (real) -my + SIM_W;	
	if (particlesCreated < MAX_PARTICLES)
	{
		real x = attractor.x + MISRENDER::frand() * 14 ;
		real y = attractor.y + MISRENDER::frand() * 14;
		if( x > r && y > r && x < SIM_W-r && y < SIM_W-r )
		{
			particles[particlesCreated].pos = glm::vec2(x,y);
			particles[particlesCreated].pos_old = glm::vec2(x,y);  
			particles[particlesCreated].force = glm::vec2(0,0);
			particles[particlesCreated].sigma = sigma;
			particles[particlesCreated].beta = beta;
			particles[particlesCreated].id = particlesCreated;

			particlesCreated++;
		}
	}
}

void FluidSolver::initParticles()
{
	// create a world with dimensions x:[-SIM_W,SIM_W] and y:[0,SIM_W*2]
	// Initialize particles
	// We will make a block of particles with a total width of 1/2 of the screen.
	particlesCreated = 0;

	real x = r;
	real y = SIM_W - r;
	for(int i=0; i<N; i++)
	{
		if(particlesCreated > N) break;

		particles[particlesCreated].pos = glm::vec2(x,y);
		particles[particlesCreated].pos_old = glm::vec2(x,y);
		particles[particlesCreated].force = glm::vec2(0,0);
		particles[particlesCreated].sigma = sigma;
		particles[particlesCreated].beta = beta;
		particles[particlesCreated].id = particlesCreated;

		if( x > SIM_W/2 )
		{
			particles[particlesCreated].color = 0;
		}else{
			particles[particlesCreated].color = 1;
		}
		

		particlesCreated++;		

		x += r * (real)0.5;

		if( x > SIM_W-r)
		{ 
			x = r;
			y -= r * 0.5f;
		}




	}
	verts.resize( N );
}

void FluidSolver::generateColorLUT()
{
	// Generate color look-up table
	double hue = 0;
	double hueInc = 6.0/255.0;
	double v = 1.0; // value 
	double s = 1.0; // saturation
	for(int i=0; i<175; i++)
	{
		double c = v * s;
		double x = c * 1.0-fabs(fmod(hue,2.0)-1.0);
		//x*=255;
		//c*=255;
#define RGB(i, r,g,b) roygbiv[i][0] = (float)r; roygbiv[i][1] = (float)g; roygbiv[i][2] = (float)b;
		if(hue >= 0 && hue <= 1){ RGB(i, c,x,0); }// red->orange
		if(hue >= 1 && hue <= 2){ RGB(i, x,c,0); }// orange->yellow
		if(hue >= 2 && hue <= 3){ RGB(i, 0,c,x); }// yellow->green
		if(hue >= 3 && hue <= 4){ RGB(i, 0,x,c); }// green->blue
		if(hue >= 4 && hue <= 5){ RGB(i, x,0,c); }// blue->indigo
		if(hue >= 5 && hue <= 6){ RGB(i, c,0,x); }// indigo->violet
#undef RGB
		hue += hueInc;
	}
}

FluidSolver::~FluidSolver()
{

	glDeleteBuffers(1, &VertexVBOID );
}

void FluidSolver::reinstallShader()
{
	//colorShader->install();
}

